% loading data
clear
tic
coarseness=150; %set 150 % parameter for grouping
minprct=10; %percentile for picking d
maxprct=90; %percentile for picking k
code0;
loaddataandgroup;
%% some other calculations from the data
ordering_list=ordering(grpsize,groupk);
fac_updatek=fac_updatek(ia2);
nfac=length(fac);
endowmenttable=table(period);
endowment=zeros(nperiod,nfac);
for j=1:nfac
    facid=fac(j);
    parfor t=1:nperiod
        periodid=period(t);
        endowment(t,j)=datafac.r(datafac.fac==facid & datafac.period==periodid);
    end
    endowmenttab_newcol=table(endowment(:,j),'VariableNames',string(facid));
    endowmenttable=[endowmenttable endowmenttab_newcol];
end
average_endowment_fac=mean(endowment,1)';

%adjusted endowment 
endowment_adj=min(max(endowment,mini_demands'),fac_updatek'+mini_demands');
average_endowment_fac_adj=mean(endowment_adj,1)';
%% get the names
Name=[biannualdata.fac_id biannualdata.name;categorical(0) 'Artifical Fac'];
Name=unique(Name,'rows');
Namee=Name(:,2);
name_agent=Namee(ismember(Name(:,1),fac));
%% emission table
emissiontable=table(period);
emission=zeros(nperiod,nfac);
for j=1:nfac
    facid=fac(j);
    for t=1:nperiod
        periodid=period(t);
        emission(t,j)=datafac.q(datafac.fac==facid & datafac.period==periodid)+mini_demands(j);
    end
    emissiontab_newcol=table(emission(:,j),'VariableNames',string(facid));
    emissiontable=[emissiontable emissiontab_newcol];
end

emp_vol=sum(abs(emission-endowment),2)/2;%-sum(abs(emission(:,57)-endowment(:,57)),2);

%%  Summary of data (for a latex table)
summarytable=zeros(56,9);
for facind=1:56
    summarytable(facind,1)=double(string(fac(facind)));
    summarytable(facind,2)=max(datafac.alloc(datafac.fac==fac(facind)));
    summarytable(facind,3)=min(datafac.alloc(datafac.fac==fac(facind)));
    summarytable(facind,4)=mean(datafac.alloc(datafac.fac==fac(facind)));
    summarytable(facind,5)=std(datafac.alloc(datafac.fac==fac(facind)));
    summarytable(facind,6)=max(datafac.emi(datafac.fac==fac(facind)));
    summarytable(facind,7)=min(datafac.emi(datafac.fac==fac(facind)));
    summarytable(facind,8)=mean(datafac.emi(datafac.fac==fac(facind)));
    summarytable(facind,9)=std(datafac.emi(datafac.fac==fac(facind)));
end
summarytable(:,2:9)=summarytable(:,2:9)./1000;
testtable=array2table(summarytable,"VariableNames",{'id','max alloc','min alloc','mean alloc','std alloc','max emission','min emission','mean emission','std emission'});
summarytable=[table(name_agent(1:56)) testtable];
%% Figure 1 plot of one firm's alloc vs emission id 12372, Mission Clay Product
figure;
bar(datafac.emi(datafac.fac=='12372'),0.7,'EdgeColor','k','FaceColor','w')
hold on;
bar(datafac.alloc(datafac.fac=='12372'),0.2,'EdgeColor','k','FaceColor','k')
legend('Initial allocations','Actual emissions')
xticks(1:nperiod)
xticklabels(period)
xtickangle(90)
xlabel('Period')
ylabel('Tons of NOx')
hold off;

%% calibration with different distribution assumptions 
warning('off','all')
powerinv;
beta_dist;
power_dist;
%%  calibration result summary 
table_noname=table(fac_updatek,average_endowment_fac,alp_fac_pi,r_fac_pi',r_fac_po',r_fac_beta');
table_noname=[table(name_agent) table_noname]
table_noname.Properties.VariableNames = {'Name','max demand','average endowment','power distribution a_i','power distribution r^*','r^* for F=1-(1-x)^a','beta r^*'};
%table_noname=table2array(table_noname);
%% Figure 2 trading volume
figure;
plot(1:nperiod,emp_vol,'LineStyle','-','Color','k','LineWidth',2.2);
hold on;
plot(1:nperiod,mean_eff_trading_vol,'LineStyle',':','Color','k','LineWidth',2.2)
plot(1:nperiod,eff_trading_vol_5,'-x','LineWidth',2,'Color','k','MarkerSize',14)
plot(1:nperiod,eff_trading_vol_95,'-o','LineWidth',2,'Color','k','MarkerSize',14)
xticks(1:10)
xticklabels(string(period))
xlabel('Period')
ylabel('Trading volume')
ax = gca;
ax.FontSize = 16; 
hold off;
%% Fig 3 revenue plots
figure;
plot(1:nperiod,emp_rev,'LineStyle','none','Marker','diamond','MarkerFaceColor','k','MarkerSize',13);
hold on;
plot(1:nperiod,r_ast_rev*ones(size(1:nperiod)),'LineStyle',':','Color','k','LineWidth',2.5,'MarkerFaceColor','w','MarkerSize',13)
plot(1:nperiod,r_sw_rev*ones(size(1:nperiod)),'-square','Color','k','LineWidth',2.5,'MarkerFaceColor','w','MarkerSize',13)
plot(1:nperiod,corner_rev*ones(size(1:nperiod)),'-o','Color','k','LineWidth',2.5,'MarkerFaceColor','w','MarkerSize',13)
xticks(1:10)
xticklabels(string(period))
xlabel('Period')
ylabel('Revenue')
ax = gca;
ax.FontSize = 16; 
hold off;
%% Fig 4 comparing endowments
figure;
bar(ea_fac_pi+mini_demands,0.8,'FaceColor',[1 1 1])
hold on
bar((r_fac_pi'+mini_demands),0.25,'FaceColor',[0 0 0])
plot(1:nfac,average_endowment_fac,'LineStyle',':','LineWidth',2.5,'Color','k');
xticks(1:nfac)
xticklabels(name_agent')
xtickangle(90)
hold off

%% K-L distances
kl_uniform=sum((average_endowment_fac./norm2).*log((norm2/57)./(average_endowment_fac)));
ea_fac_pi_nonzero=ea_fac_pi+mini_demands;
ea_fac_pi_nonzero(ea_fac_pi+mini_demands==0)=average_endowment_fac(ea_fac_pi+mini_demands==0);
kl_ea=sum((average_endowment_fac./norm2).*log((ea_fac_pi_nonzero)./(average_endowment_fac)));
r_fac_pi_nonzero=r_fac_pi'+mini_demands;
r_fac_pi_nonzero(r_fac_pi'+mini_demands==0)=average_endowment_fac(r_fac_pi'+mini_demands==0);
kl_rast=sum((average_endowment_fac./norm2).*log((r_fac_pi_nonzero)./(average_endowment_fac)));

%% K-L divergence for assuming different distributions %15 is the zero entry
r_fac_po_nonzero=r_fac_po';
zeroentry=find(~r_fac_po_nonzero)
r_fac_po_nonzero(zeroentry)=[];
r_fac_pi_nonzero=r_fac_pi';
r_fac_pi_nonzero(zeroentry)=[];
r_fac_beta_nonzero=r_fac_beta';
r_fac_beta_nonzero(zeroentry)=[];
KL_po=sum(r_fac_pi_nonzero.*log((r_fac_po_nonzero)./(r_fac_pi_nonzero)))
KL_beta=sum(r_fac_pi_nonzero.*log((r_fac_beta_nonzero)./(r_fac_pi_nonzero)))

%% Fig OA. 
figure;scatter(fac_indivk,mean(emission),'filled');hl=lsline;set(hl(1),'color','r');xlabel('Maximum emission');hold off
figure;
histogram(fac_indivk,40);
xlabel('Maximum emission');
toc

